%% integral_loc
function [points, weights] = integral_loc(V, int_type)
[int_points_ref, int_weights_ref] = integral_ref(int_type);
if 2 == size(V,1) && strcmp(int_type, "1d")
    AB = 0.5*(V(2,:)-V(1,:));
    weights = norm(AB)*int_weights_ref;
    points = AB.*(int_points_ref+1)+V(1,:);
elseif 3 == size(V,1) && strcmp(int_type, "2d")
    Jacobi = V(2:3,:)-V(1,:);
    J_det = det(Jacobi);
    weights = abs(J_det)*int_weights_ref;
    points = int_points_ref*Jacobi+V(1,:);
else
    error("Invalid gauss-lengdre type.");
end
end

%% integral_ref
function [points, weights] = integral_ref(int_type)
switch int_type
    case "1d" % gauss-lengdre(3,4,9)
        % points = [-0.7745966692414833770358531;
        %     0;
        %     0.7745966692414833770358531];
        % weights = [0.5555555555555555555555556, 0.8888888888888888888888889, 0.5555555555555555555555556];

        points = [-0.8611363115940525752239465;
            -0.3399810435848562648026658;
            0.3399810435848562648026658;
            0.8611363115940525752239465];
        weights = [0.3478548451374538573730639, 0.6521451548625461426269361, 0.6521451548625461426269361, 0.3478548451374538573730639];
        
        % points = [-0.9681602395076260898355762;
        %     -0.8360311073266357942994298;
        %     -0.6133714327005903973087020;
        %     -0.3242534234038089290385380;
        %     0;
        %     0.3242534234038089290385380;
        %     0.6133714327005903973087020;
        %     0.8360311073266357942994298;
        %     0.9681602395076260898355762];
        % weights = [0.0812743883615744119718922, 0.180648160694

    case "2d" % Dunavant(7)
        points = [1/3, 1/3;
            0.470142064105115, 0.470142064105115;
            0.470142064105115, 0.059715871789770;
            0.059715871789770, 0.470142064105115;
            0.101286507323456, 0.101286507323456;
            0.101286507323456, 0.797426985353087;
            0.797426985353087, 0.101286507323456];
        weights = [0.1125, 0.066197076394253, 0.066197076394253, 0.066197076394253, 0.062969590272414, 0.062969590272414, 0.062969590272414];

        % points = [0.11270166537925829786104259255808, 0.1;
        %     0.11270166537925829786104259255808, 0.44364916731037085106947870372096;
        %     0.11270166537925829786104259255808, 0.78729833462074172434341789994505;
        %     0.5, 0.056350832689629148930521296279039;
        %     0.5, 0.25;
        %     0.5, 0.44364916731037085106947870372096;
        %     0.88729833462074170213895740744192, 0.01270166537925830792243875322356;
        %     0.88729833462074170213895740744192, 0.056350832689629148930521296279039;
        %     0.88729833462074170213895740744192, 0.1];
        % weights = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];

        % pts_1d = [-0.861136311594053;
        %     -0.339981043584856;
        %     0.339981043584856;
        %     0.861136311594053];
        % wts_1d = [0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454];
        % [xi_rect, eta_rect] = meshgrid(pts_1d, pts_1d);
        % xi_rect = xi_rect(:);
        % eta_rect = eta_rect(:);
        % [w_xi, w_eta] = meshgrid(wts_1d, wts_1d);
        % weights_rect = w_xi(:) .* w_eta(:);
        % xi_tri  = 0.5*(1 + xi_rect);
        % eta_tri = 0.25*(1 - xi_rect).*(1 + eta_rect);
        % J = 1/8*(1 - xi_rect);
        % points = [xi_tri, eta_tri];
        % weights = (weights_rect.*J)';

    otherwise
        error("Invalid gauss-lengdre type.");
end
end